{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3075a55c",
   "metadata": {},
   "source": [
    "@authors: dveske and zmarka"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1709957e",
   "metadata": {},
   "source": [
    "Joint skymaps in paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5f2356a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import json\n",
    "import numpy\n",
    "from matplotlib import pyplot as plt\n",
    "from hpmoc import PartialUniqSkymap\n",
    "from hpmoc.plot import get_wcs, plot, gridplot\n",
    "from hpmoc.points import PointsTuple\n",
    "\n",
    "FITS = '1262142545.615/lvc_skymap.fits.gz' #CBC event in Figure 1\n",
    "NUS = '1262142545.615/icecube_neutrino_list.json'\n",
    "\n",
    "# FITS = '1241247887.938/lvc_skymap.fits.gz' #cWB event in Figure 3\n",
    "# NUS = '1241247887.938/icecube_neutrino_list.json'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac5fb92",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(NUS) as f:\n",
    "    d = json.load(f)\n",
    "\n",
    "neutrinos = []\n",
    "\n",
    "for n in d:\n",
    "    coords = (n['ra'], n['dec'], n['sigma'])\n",
    "    neutrinos.append(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed827c2c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "points = PointsTuple(neutrinos)\n",
    "\n",
    "m = PartialUniqSkymap.read(FITS,strategy='ligo')\n",
    "\n",
    "m.plot(points, fig={'dpi': 100, 'figsize': (8, 4)},\n",
    "                             projection='MOL', sigmas=[],\n",
    "                             cbar={'orientation': 'vertical',\n",
    "                                   'shrink': 0.8, 'pad': 0.06},\n",
    "                             scatter_marker_size=60, scatter_label_size=12)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e09edac",
   "metadata": {},
   "source": [
    "Figure 5. in paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6210df42",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=numpy.genfromtxt('gwhen_o3_UL.csv',delimiter=',')\n",
    "\n",
    "count=0\n",
    "plt.figure()\n",
    "for typ in['BNS','NSBH','BBH']:\n",
    "    if(typ=='BNS'):\n",
    "        color='tab:blue'\n",
    "    elif(typ=='NSBH'):\n",
    "        color='tab:orange'\n",
    "    else:\n",
    "        color='tab:green'\n",
    "    for fb in [1,100]:\n",
    "        ul=data[count]  \n",
    "        if(fb==1):\n",
    "            plt.plot(numpy.logspace(47,55,100),ul,color=color,label=typ)\n",
    "        elif(fb==100):\n",
    "            plt.plot(numpy.logspace(47,55,100),ul,'--',color=color)\n",
    "        count+=1\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.ylabel('90% upper limit for GWHEN events [Gpc$^{-3}$ yr$^{-1}$]')\n",
    "plt.xlabel(r'Total $E_{\\nu}$ [erg]')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.fill_between(numpy.logspace(47,55,100),10,1700,color='none',edgecolor='tab:blue',linestyle='dotted',alpha=1,linewidth=2)\n",
    "plt.fill_between(numpy.logspace(47,55,100),7.8,140,color='none',edgecolor='tab:orange',linestyle='dotted',alpha=1,linewidth=2)\n",
    "plt.fill_between(numpy.logspace(47,55,100),17.9,44,color='none',edgecolor='tab:green',linestyle='dotted',alpha=1,linewidth=2)\n",
    "plt.text(10**48, 500, 'O3 BNS',color='tab:blue')\n",
    "plt.text(10**49.5, 65, 'O3 NSBH',color='tab:orange')\n",
    "plt.text(10**51, 20, 'O3 BBH',color='tab:green')\n",
    "plt.title(r'Solid: $f_{\\rm beam}$=1, Dashed: $f_{\\rm beam}$=100')\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
